C     ALGORITHM 476 COLLECTED ALGORITHMS FROM ACM.                              
C     ALGORITHM APPEARED IN COMM. ACM, VOL. 17, NO. 04,                         
C     P. 220.                                                                   
      SUBROUTINE CURV1(N, X, Y, SLP1, SLPN, YP, TEMP, SIGMA)                    
      INTEGER N                                                                 
      REAL X(N), Y(N), SLP1, SLPN, YP(N), TEMP(N), SIGMA                        
C THIS SUBROUTINE DETERMINES THE PARAMETERS NECESSARY TO                        
C COMPUTE AN INTERPOLATORY SPLINE UNDER TENSION THROUGH                         
C A SEQUENCE OF FUNCTIONAL VALUES. THE SLOPES AT THE TWO                        
C ENDS OF THE CURVE MAY BE SPECIFIED OR OMITTED. FOR ACTUAL                     
C COMPUTATION OF POINTS ON THE CURVE IT IS NECESSARY TO CALL                    
C THE FUNCTION CURV2.                                                           
C ON INPUT--                                                                    
C N IS THE NUMBER OF VALUES TO BE INTERPOLATED (N.GE.2),                        
C X IS AN ARRAY OF THE N INCREASING ABSCISSAE OF THE                            
C FUNCTIONAL VALUES,                                                            
C Y IS AN ARRAY OF THE N ORDINATES OF THE VALUES,(I.E.Y(K)                      
C IS THE FUNCTIONAL VALUE CORRESPONDING TO X(K)),                               
C SLP1 AND SLPN CONTAIN THE DESIRED VALUES FOR THE FIRST                        
C DERIVATIVE OF THE CURVE AT X(1) AND X(N), RESPECTIVELY.                       
C IF THE QUANTITY SIGMA IS NEGATIVE THESE VALUES WILL BE                        
C DETERMINED INTERNALLY AND THE USER NEED ONLY FURNISH                          
C PLACE-HOLDING PARAMETERS FOR SLP1 AND SLPN. SUCH PLACE-                       
C HOLDING PARAMETERS WILL BE IGNORED BUT NOT DESTROYED,                         
C YP IS AN ARRAY OF LENGTH AT LEAST N                                           
C TEMP IS AN ARRAY OF LENGTH AT LEAST N WHICH IS USED FOR                       
C SCRATCH STORAGE,                                                              
C AND                                                                           
C SIGMA CONTAINS THE TENSION FACTOR. THIS IS NON-ZERO AND                       
C INDICATES THE CURVINESS DESIRED. IF ABS(SIGMA) IS NEARLY                      
C ZERO (E.G. .001) THE RESULTING CURVE IS APPROXIMATELY A                       
C CUBIC SPLINE. IF ABS(SIGMA) IS LARGE (E.G. 50.) THE                           
C RESULTING CURVE IS NEARLY A POLYGONAL LINE. THE SIGN                          
C OF SIGMA INDICATES WHETHER THE DERIVATIVE INFORMATION                         
C HAS BEEN INPUT OR NOT. IF SIGMA IS NEGATIVE THE ENDPOINT                      
C DERIVATIVES WILL BE DETERMINED INTERNALLY. A STANDARD                         
C VALUE FOR SIGMA IS APPROXIMATELY 1. IN ABSOLUTE VALUE.                        
C ON OUTPUT--                                                                   
C YP CONTAINS VALUES PROPORTIONAL TO THE SECOND DERIVATIVE                      
C OF THE CURVE AT THE GIVEN NODES.                                              
C N,X,Y,SLP1,SLPN AND SIGMA ARE UNALTERED,                                      
      NM1 = N - 1                                                               
      NP1 = N + 1                                                               
      DELX1 = X(2) - X(1)                                                       
      DX1 = (Y(2)-Y(1))/DELX1                                                   
C DETERMINE SLOPES IF NECESSARY                                                 
      IF (SIGMA.LT.0.) GO TO 50                                                 
      SLPP1 = SLP1                                                              
      SLPPN = SLPN                                                              
C DENORMALIZE TENSION FACTOR                                                    
   10 SIGMAP = ABS(SIGMA)*FLOAT(N-1)/(X(N)-X(1))                                
C SET UP RIGHT HAND SIDE AND TRIDIAGONAL SYSTEM FOR YP AND                      
C PERFORM FORWARD ELIMINATION                                                   
      DELS = SIGMAP*DELX1                                                       
      EXPS = EXP(DELS)                                                          
      SINHS = .5*(EXPS-1./EXPS)                                                 
      SINHIN = 1./(DELX1*SINHS)                                                 
      DIAG1 = SINHIN*(DELS*.5*(EXPS+1./EXPS)-SINHS)                             
      DIAGIN = 1./DIAG1                                                         
      YP(1) = DIAGIN*(DX1-SLPP1)                                                
      SPDIAG = SINHIN*(SINHS-DELS)                                              
      TEMP(1) = DIAGIN*SPDIAG                                                   
      IF (N.EQ.2) GO TO 30                                                      
      DO 20 I=2,NM1                                                             
        DELX2 = X(I+1) - X(I)                                                   
        DX2 = (Y(I+1)-Y(I))/DELX2                                               
        DELS = SIGMAP*DELX2                                                     
        EXPS = EXP(DELS)                                                        
        SINHS = .5*(EXPS-1./EXPS)                                               
        SINHIN = 1./(DELX2*SINHS)                                               
        DIAG2 = SINHIN*(DELS*(.5*(EXPS+1./EXPS))-SINHS)                         
        DIAGIN = 1./(DIAG1+DIAG2-SPDIAG*TEMP(I-1))                              
        YP(I) = DIAGIN*(DX2-DX1-SPDIAG*YP(I-1))                                 
        SPDIAG = SINHIN*(SINHS-DELS)                                            
        TEMP(I) = DIAGIN*SPDIAG                                                 
        DX1 = DX2                                                               
        DIAG1 = DIAG2                                                           
   20 CONTINUE                                                                  
   30 DIAGIN = 1./(DIAG1-SPDIAG*TEMP(NM1))                                      
      YP(N) = DIAGIN*(SLPPN-DX2-SPDIAG*YP(NM1))                                 
C PERFORM BACK SUBSTITUTION                                                     
      DO 40 I=2,N                                                               
        IBAK = NP1 - I                                                          
        YP(IBAK) = YP(IBAK) - TEMP(IBAK)*YP(IBAK+1)                             
   40 CONTINUE                                                                  
      RETURN                                                                    
   50 IF (N.EQ.2) GO TO 60                                                      
C IF NO DERIVATIVES ARE GIVEN USE SECOND ORDER POLYNOMIAL                       
C INTERPOLATION ON INPUT DATA FOR VALUES AT ENDPOINTS.                          
      DELX2 = X(3) - X(2)                                                       
      DELX12 = X(3) - X(1)                                                      
      C1 = -C(DELX12+DELX1)/DELX12/DELX1                                        
      C2 = DELX12/DELX1/DELX2                                                   
      C3 = -DELX1/DELX12/DELX2                                                  
      SLPP1 = C1*Y(1) + C2*Y(2) + C3*Y(3)                                       
      DELN = X(N) - X(NM1)                                                      
      DELNM1 = X(NM1) - X(N-2)                                                  
      DELNN = X(N) - X(N-2)                                                     
      C1 = (DELNN+DELN)/DELNN/DELN                                              
      C2 = -DELNN/DELN/DELNM1                                                   
      C3 = DELN/DELNN/DELNM1                                                    
      SLPPN = C3*Y(N-2) + C2*Y(NM1) + C1*Y(N)                                   
      GO TO 10                                                                  
C IF ONLY TWO POINTS AND NO DERIVATIVES ARE GIVEN, USE                          
C STRAIGHT LINE FOR CURVE                                                       
   60 YP(1) = 0.                                                                
      YP(2) = 0.                                                                
      RETURN                                                                    
      END                                                                       
      FUNCTION CURV2(T, N, X, Y, YP, SIGMA, IT)                                 
      INTEGER N, IT                                                             
      REAL T, X(N), Y(N), YP(N), SIGMA                                          
C THIS FUNCTION INTERPOLATES A CURVE AT A GIVEN POINT                           
C USING A SPLINE UNDER TENSION. THE SUBROUTINE CURV1 SHOULD                     
C BE CALLED EARLIER TO DETERMINE CERTAIN NECESSARY                              
C PARAMETERS.                                                                   
C ON INPUT--                                                                    
C T CONTAINS A REAL VALUE TO BE MAPPED ONTO THE INTERPO-                        
C LATING CURVE.                                                                 
C N CONTAINS THE NUMBER OF POINTS WHICH WERE INTERPOLATED                       
C TO DETERMINE THE CURVE,                                                       
C X AND Y ARE ARRAYS CONTAINING THE ORDINATES AND ABCISSAS                      
C OF THE INTERPOLATED POINTS,                                                   
C YP IS AN ARRAY WITH VALUES PROPORTIONAL TO THE SECOND                         
C DERIVATIVE OF THE CURVE AT THE NODES                                          
C SIGMA CONTAINS THE TENSION FACTOR (ITS SIGN IS IGNORED)                       
C IT IS AN INTEGER SWITCH. IF IT IS NOT 1 THIS INDICATES                        
C THAT THE FUNCTION HAS BEEN CALLED PREVIOUSLY (WITH N,X,                       
C Y,YP, AND SIGMA UNALTERED) AND THAT THIS VALUE OF T                           
C EXCEEDS THE PREVIOUS VALUE. WITH SUCH INFORMATION THE                         
C FUNCTION IS ABLE TO PERFORM THE INTERPOLATION MUCH MORE                       
C RAPIDLY. IF A USER SEEKS TO INTERPOLATE AT A SEQUENCE                         
C OF POINTS, EFFICIENCY IS GAINED BY ORDERING THE VALUES                        
C INCREASING AND SETTING IT TO THE INDEX OF THE CALL.                           
C IF IT IS 1 THE SEARCH FOR THE INTERVAL (X(K),X(K+1))                          
C CONTAINING T STARTS WITH K=1.                                                 
C THE PARAMETERS N,X,Y,YP AND SIGMA SHOULD BE INPUT                             
C UNALTERED FROM THE OUTPUT OF CURV1.                                           
C ON OUTPUT--                                                                   
C CURV2 CONTAINS THE INTERPOLATED VALUE. FOR T LESS THAN                        
C X(1) CURV2 = Y(1). FOR T GREATER THAN X(N) CURV2 = Y(N).                      
C NONE OF THE INPUT PARAMETERS ARE ALTERED.                                     
      S = X(N) - X(1)                                                           
C DENORMALIZE SIGMA                                                             
      SIGMAP = ABS(SIGMA)*FLOAT(N-1)/S                                          
C IF IT.NE.1 START SEARCH WHERE PREVIOUSLY TERMINATED,                          
C OTHERWISE START FROM BEGINNING                                                
      IF (IT.EQ.1) I1 = 2                                                       
C SEARCH FOR INTERVAL                                                           
   10 DO 20 I=I1,N                                                              
        IF (X(I)-T) 20, 20, 30                                                  
   20 CONTINUE                                                                  
      I = N                                                                     
C CHECK TO INSURE CORRECT INTERVAL                                              
   30 IF (X(I-1).LE.T .OR. T.LE.X(1)) GO TO 40                                  
C RESTART SEARCH  AND RESET I1                                                  
C ( INPUT ""IT"" WAS INCORRECT )                                                
      I1 = 2                                                                    
      GO TO 10                                                                  
C SET UP AND PERFORM INTERPOLATION                                              
   40 DEL1 = T - X(I-1)                                                         
      DEL2 = X(I) - T                                                           
      DELS = X(I) - X(I-1)                                                      
      EXPS1 = EXP(SIGMAP*DEL1)                                                  
      SINHD1 = .5*(EXPS1-1./EXPS1)                                              
      EXPS = EXP(SIGMAP*DEL2)                                                   
      SINHD2 = .5*(EXPS-1./EXPS)                                                
      EXPS = EXPS1*EXPS                                                         
      SINHS = .5*(EXPS-1./EXPS)                                                 
      CURV2 = (YP(I)*SINHD1+YP(I-1)*SINHD2)/SINHS +                             
     * ((Y(I)-YP(I))*DEL1+(Y(I-1)-YP(I-1))*DEL2)/DELS                           
      I1 = I                                                                    
      RETURN                                                                    
      END                                                                       
      SUBROUTINE KURV1(N, X, Y, SLP1, SLPN, XP, YP, TEMP, S,                    
     * SIGMA)                                                                   
C THIS SUBROUTINE DETERMINES THE PARAMETERS NECESSARY TO                        
C COMPUTE A SPLINE UNDER TENSION PASSING THROUGH A SEQUENCE                     
C OF PAIRS  (X(1),Y(1)),...,(X(N),Y(N)) IN THE PLANE. THE                       
C SLOPES AT THE TWO ENDS OF THE CURVE MAY BE SPECIFIED OR                       
C OMITTED. FOR ACTUAL COMPUTATION OF POINTS ON THE CURVE IT                     
C IS NECESSARY TO CALL THE SUBROUTINE KURV2.                                    
C ON INPUT--                                                                    
C N IS THE NUMBER OF POINTS TO BE INTERPOLATED (N.GE.2),                        
C X IS AN ARRAY CONTAINING THE N X-COORDINATES OF THE                           
C POINTS,                                                                       
C Y IS AN ARRAY CONTAINING THE N Y-COORDINATES OF THE                           
C POINTS,                                                                       
C SLP1 AND SLPN CONTAIN THE DESIRED VALUES FOR THE SLOPE                        
C OF THE CURVE AT  (X(1),Y(1)) AND  (X(N),Y(N)), RESPEC-                        
C TIVELY. THESE QUANTITIES ARE IN DEGREES AND MEASURED                          
C COUNTERCLOCKWISE FROM THE POSITIVE X-AXIS. THE POSITIVE                       
C SENSE OF THE CURVE IS ASSUMED TO BE THAT MOVING FROM THE                      
C POINT 1 TO POINT N. IF THE QUANTITY SIGMA IS NEGATIVE                         
C THESE SLOPES WILL BE DETERMINED INTERNALLY AND THE USER                       
C NEED ONLY FURNISH PLACE-HOLDING PARAMETERS FOR SLP1 AND                       
C SLPN. SUCH PLACE-HOLDING PARAMETERS WILL BE IGNORED BUT                       
C NOT DESTROYED,                                                                
C XP,YP ARE ARRAYS OF LENGTH AT LEAST N,                                        
C TEMP IS AN ARRAY OF LENGTH AT LEAST N WHICH IS USED FOR                       
C SCRATCH STORAGE,                                                              
C AND                                                                           
C SIGMA CONTAINS THE TENSION FACTOR. THIS IS NON-ZERO AND                       
C INDICATES THE CURVINESS DESIRED. IF ABS(SIGMA) IS VERY                        
C LARGE (E.G. 50.) THE RESULTING CURVE IS VERY NEARLY A                         
C POLYGONAL LINE. THE SIGN OF SIGMA INDICATES WHETHER                           
C SLOPE IN FORMATION HAS BEEN INPUT OR NOT. IF SIGMA IS                         
C NEGATIVE THE END-POINT SLOPES WILL BE DETERMINED                              
C INTERNALLY. A STANDARD VALUE FOR SIGMA IS APPROXIMATELY                       
C 1. IN ABSOLUTE VALUE.                                                         
C ON OUTPUT--                                                                   
C N,X,Y,SLP1,SLPN, AND SIGMA ARE UNALTERED,                                     
C XP AND YP CONTAIN INFORMATION ABOUT THE CURVATURE OF THE                      
C CURVE AT THE GIVEN NODES,                                                     
C AND                                                                           
C S CONTAINS THE POLYGONAL ARCLENGTH OF THE CURVE.                              
      INTEGER N                                                                 
      REAL X(N), Y(N), XP(N), YP(N), TEMP(N), S, SIGMA                          
      DEGRAD = 3.1415926535897932/180.                                          
      NM1 = N - 1                                                               
      NP1 = N + 1                                                               
      DELX1 = X(2) - X(1)                                                       
      DELY1 = Y(2) - Y(1)                                                       
      DELS1 = SQRT(DELX1*DELX1+DELY1*DELY1)                                     
      DX1 = DELX1/DELS1                                                         
      DY1 = DELY1/DELS1                                                         
C DETERMINE SLOPES IF NECESSARY                                                 
      IF (SIGMA.LT.0.) GO TO 70                                                 
      SLPP1 = SLP1*DEGRAD                                                       
      SLPPN = SLPN*DEGRAD                                                       
C SET UP RIGHT HAND SIDES OF TRIDIAGONAL LINEAR SYSTEM FOR XP                   
C AND YP                                                                        
   10 XP(1) = DX1 - COS(SLPP1)                                                  
      YP(1) = DY1 - SIN(SLPP1)                                                  
      TEMP(1) = DELS1                                                           
      S = DELS1                                                                 
      IF (N.EQ.2) GO TO 30                                                      
      DO 20 I=2,NM1                                                             
        DELX2 = X(I+1) - X(I)                                                   
        DELY2 = Y(I+1) - Y(I)                                                   
        DELS2 = SQRT(DELX2*DELX2+DELY2*DELY2)                                   
        DX2 = DELX2/DELS2                                                       
        DY2 = DELY2/DELS2                                                       
        XP(I) = DX2 - DX1                                                       
        YP(I) = DY2 - DY1                                                       
        TEMP(I) = DELS2                                                         
        DELX1 = DELX2                                                           
        DELY1 = DELY2                                                           
        DELS1 = DELS2                                                           
        DX1 = DX2                                                               
        DY1 = DY2                                                               
C ACCUMULATE POLYGONAL ARCLENGTH                                                
        S = S + DELS1                                                           
   20 CONTINUE                                                                  
   30 XP(N) = COS(SLPPN) - DX1                                                  
      YP(N) = SIN(SLPPN) - DY1                                                  
C DENORMALIZE TENSION FACTOR                                                    
      SIGMAP = ABS(SIGMA)*FLOAT(N-1)/S                                          
C PERFORM FORWARD ELIMINATION ON TRIDIAGONAL SYSTEM                             
      DELS = SIGMAP*TEMP(1)                                                     
      EXPS = EXP(DELS)                                                          
      SINHS = .5*(EXPS-1./EXPS)                                                 
      SINHIN = 1./(TEMP(1)*SINHS)                                               
      DIAG1 = SINHIN*(DELS*.5*(EXPS+1./EXPS)-SINHS)                             
      DIAGIN = 1./DIAG1                                                         
      XP(1) = DIAGIN*XP(1)                                                      
      YP(1) = DIAGIN*YP(1)                                                      
      SPDIAG = SINHIN*(SINHS-DELS)                                              
      TEMP(1) = DIAGIN*SPDIAG                                                   
      IF (N.EQ.2) GO TO 50                                                      
      DO 40 I=2,NM1                                                             
        DELS = SIGMAP*TEMP(I)                                                   
        EXPS = EXP(DELS)                                                        
        SINHS = .5*(EXPS-1./EXPS)                                               
        SINHIN = 1./(TEMP(I)*SINHS)                                             
        DIAG2 = SINHIN*(DELS*(.5*(EXPS+1./EXPS))-SINHS)                         
        DIAGIN = 1./(DIAG1+DIAG2-SPDIAG*TEMP(I-1))                              
        XP(I) = DIAGIN*(XP(I)-SPDIAG*XP(I-1))                                   
        YP(I) = DIAGIN*(YP(I)-SPDIAG*YP(I-1))                                   
        SPDIAG = SINHIN*(SINHS-DELS)                                            
        TEMP(I) = DIAGIN*SPDIAG                                                 
        DIAG1 = DIAG2                                                           
   40 CONTINUE                                                                  
   50 DIAGIN = 1./(DIAG1-SPDIAG*TEMP(NM1))                                      
      XP(N) = DIAGIN*(XP(N)-SPDIAG*XP(NM1))                                     
      YP(N) = DIAGIN*(YP(N)-SPDIAG*YP(NM1))                                     
C PERFORM BACK SUBSTITUTION                                                     
      DO 60 I=2,N                                                               
        IBAK = NP1 - I                                                          
        XP(IBAK) = XP(IBAK) - TEMP(IBAK)*XP(IBAK+1)                             
        YP(IBAK) = YP(IBAK) - TEMP(IBAK)*YP(IBAK+1)                             
   60 CONTINUE                                                                  
      RETURN                                                                    
   70 IF (N.EQ.2) GO TO 80                                                      
C IF NO SLOPES ARE GIVEN, USE SECOND ORDER INTERPOLATION ON                     
C INPUT DATA FOR SLOPES AT ENDPOINTS                                            
      DELS2 = SQRT((X(3)-X(2))**2+(Y(3)-Y(2))**2)                               
      DELS12 = DELS1 + DELS2                                                    
      C1 = -C(DELS12+DELS1)/DELS12/DELS1                                        
      C2 = DELS12/DELS1/DELS2                                                   
      C3 = -DELS1/DELS12/DELS2                                                  
      SX = C1*X(1) + C2*X(2) + C3*X(3)                                          
      SY = C1*Y(1) + C2*Y(2) + C3*Y(3)                                          
      SLPP1 = ATAN2(SY,SX)                                                      
      DELNM1 = SQRT((X(N-2)-X(NM1))**2+(Y(N-2)-Y(NM1))**2)                      
      DELN = SQRT((X(NM1)-X(N))**2+(Y(NM1)-Y(N))**2)                            
      DELNN = DELNM1 + DELN                                                     
      C1 = (DELNN+DELN)/DELNN/DELN                                              
      C2 = -DELNN/DELN/DELNM1                                                   
      C3 = DELN/DELNN/DELNM1                                                    
      SX = C3*X(N-2) + C2*X(NM1) + C1*X(N)                                      
      SY = C3*Y(N-2) + C2*Y(NM1) + C1*Y(N)                                      
      SLPPN = ATAN2(SY,SX)                                                      
      GO TO 10                                                                  
C IF ONLY TWO POINTS AND NO SLOPES ARE GIVEN, USE STRAIGHT                      
C LINE SEGMENT FOR CURVE                                                        
   80 XP(1) = 0.                                                                
      XP(2) = 0.                                                                
      YP(1) = 0.                                                                
      YP(2) = 0.                                                                
      RETURN                                                                    
      END                                                                       
      SUBROUTINE KURV2(T, XS, YS, N, X, Y, XP, YP, S, SIGMA)                    
      INTEGER N                                                                 
      REAL T, XS, YS, X(N), Y(N), XP(N), YP(N), S, SIGMA                        
C THIS SUBROUTINE PERFORMS THE MAPPING OF POINTS IN THE                         
C INTERVAL (0.,1.) ONTO A CURVE IN THE PLANE. THE SUBROUTINE                    
C KURV1 SHOULD BE CALLED EARLIER TO DETERMINE CERTAIN                           
C NECESSARY PARAMETERS. THE RESULTING CURVE HAS A PARAMETRIC                    
C REPRESENTATION BOTH OF WHOSE COMPONENTS ARE SPLINES UNDER                     
C TENSION AND FUNCTIONS OF THE POLYGONAL ARCLENGTH PARAMETER                    
C ON INPUT--                                                                    
C T CONTAINS A REAL VALUE OF ABSOLUTE VALUE LESS THAN OR                        
C EQUAL TO 1. TO BE MAPPED TO A POINT ON THE CURVE. THE                         
C SIGN OF T IS IGNORED AND THE INTERVAL (0.,1.) IS MAPPED                       
C ONTO THE ENTIRE CURVE. IF T IS NEGATIVE THIS INDICATES                        
C THAT THE SUBROUTINE HAS BEEN CALLED PREVIOUSLY (WITH ALL                      
C OTHER INPUT VARIABLES UNALTERED) AND THAT THIS VALUE OF                       
C T EXCEEDS THE PREVIOUS VALUE IN ABSOLUTE VALUE. WITH                          
C SUCH INFORMATION THE SUBROUTINE IS ABLE TO MAP THE POINT                      
C MUCH MORE RAPIDLY. THUS IF THE USER SEEKS TO MAP A                            
C SEQUENCE OF POINTS ONTO THE SAME CURVE, EFFICIENCY IS                         
C GAINED BY ORDERING THE VALUES INCREASING IN MAGNITUDE                         
C AND SETTING THE SIGNS OF ALL BUT THE FIRST, NEGATIVE,                         
C N CONTAINS THE NUMBER OF POINTS WHICH WERE INTERPOLATED                       
C TO DETERMINE THE CURVE,                                                       
C X AND Y ARE ARRAYS CONTAINING THE X- AND Y-COORDINATES                        
C OF THE INTERPOLATED POINTS,                                                   
C XP AND YP ARE THE ARRAYS OUTPUT FROM KURV2 CONTAINING                         
C CURVATURE INFORMATION,                                                        
C S CONTAINS THE POLYGONAL ARCLENGTH OF THE CURVE,                              
C SIGMA CONTAINS THE TENSION FACTOR (ITS SIGN IS IGNORED).                      
C THE PARAMETERS N,X,Y,XP,YP,S,AND SIGMA SHOULD BE INPUT                        
C UNALTERED FROM THE OUTPUT OF KURV1.                                           
C ON OUTPUT--                                                                   
C XS AND YS CONTAIN THE X-ANDY-COORDINATES OF THE IMAGE                         
C POINT ON THE CURVE.                                                           
C T,N,X,Y,XP,YP,S, AND SIGMA ARE UNALTERED.                                     
C DENORMALIZE SIGMA                                                             
      SIGMAP = ABS(SIGMA)*FLOAT(N-1)/S                                          
C STRETCH UNIT INTERVAL INTO ARCLENGTH DISTANCE                                 
      TN = ABS(T*S)                                                             
C FOR NEGATIVE T START SEARCH WHERE PREVIOUSLY TERMINATED,                      
C OTHERWISE START FROM BEGINNING                                                
      IF (T.LT.0.) GO TO 10                                                     
      I1 = 2                                                                    
      XS = X(1)                                                                 
      YS = Y(1)                                                                 
      SUM = 0.                                                                  
      IF (T.LE.0.) RETURN                                                       
   10 CONTINUE                                                                  
C DETERMINE INTO WHICH SEGMENT TN IS MAPPED                                     
      DO 30 I=I1,N                                                              
        DELX = X(I) - X(I-1)                                                    
        DELY = Y(I) - Y(I-1)                                                    
        DELS = SQRT(DELX*DELX+DELY*DELY)                                        
        IF (SUM+DELS-TN) 20, 40, 40                                             
   20   SUM = SUM + DELS                                                        
   30 CONTINUE                                                                  
C IF ABS(T) IS GREATER THAN 1., RETURN TERMINAL POINT ON                        
C CURVE                                                                         
      XS = X(N)                                                                 
      YS = Y(N)                                                                 
      RETURN                                                                    
C SET UP AND PERFORM INTERPOLATION                                              
   40 DEL1 = TN - SUM                                                           
      DEL2 = DELS - DEL1                                                        
      EXPS1 = EXP(SIGMAP*DEL1)                                                  
      SINHD1 = .5*(EXPS1-1./EXPS1)                                              
      EXPS = EXP(SIGMAP*DEL2)                                                   
      SINHD2 = .5*(EXPS-1./EXPS)                                                
      EXPS = EXPS1*EXPS                                                         
      SINHS = .5*(EXPS-1./EXPS)                                                 
      XS = (XP(I)*SINHD1+XP(I-1)*SINHD2)/SINHS +                                
     * ((X(I)-XP(I))*DEL1+(X(I-1)-XP(I-1))*DEL2)/DELS                           
      YS = (YP(I)*SINHD1+YP(I-1)*SINHD2)/SINHS +                                
     * ((Y(I)-YP(I))*DEL1+(Y(I-1)-YP(I-1))*DEL2)/DELS                           
      I1 = I                                                                    
      RETURN                                                                    
      END                                                                       
      SUBROUTINE KURVP1(N, X, Y, XP, YP, TEMP, S, SIGMA)                        
      INTEGER N                                                                 
      REAL X(N), Y(N), XP(N), YP(N), TEMP(1), S, SIGMA                          
C THIS SUBROUTINE DETERMINES THE PARAMETERS NECESSARY TO                        
C COMPUTE A SPLINE UNDER TENSION FORMING A CLOSED CURVE IN                      
C THE PLANE AND PASSING THROUGH A SEQUENCE OF PAIRS                             
C (X(1),Y(1)),...,(X(N),Y(N)). FOR ACTUAL COMPUTATION OF                        
C POINTS ON THE CURVE IT IS NECESSARY TO CALL THE SUBROUTINE                    
C KURVP2.                                                                       
C ON INPUT--                                                                    
C N IS THE NUMBER OF POINTS TO BE INTERPOLATED (N.GE.2),                        
C X IS AN ARRAY CONTAINING THE N X-COORDINATES OF THE                           
C POINTS,                                                                       
C Y IS AN ARRAY CONTAINING THE N Y-COORDINATES OF THE                           
C POINTS,                                                                       
C XP,YP ARE ARRAYS OF LENGTH AT LEAST N,                                        
C TEMP IS AN ARRAY OF LENGTH AT LEAST 2*N WHICH IS USED                         
C FOR SCRATCH STORAGE,                                                          
C AND                                                                           
C SIGMA CONTAINS THE TENSION FACTOR. THIS IS A NON-ZERO                         
C QUANTITY (WHOSE SIGN IS IGNORED) WHICH INDICATES THE                          
C CURVINESS DESIRED. IF ABS(SIGMA) IS VERY LARGE (E.G. 50.                      
C ) THE RESULTING CURVE IS VERY A POLYGON. A STANDARD                           
C VALUE FOR SIGMA IS APPROXIMATELY 1. IN ABSOLUTE VALUE.                        
C ON OUTPUT--                                                                   
C N,X,Y, AND SIGMA ARE UNALTERED,                                               
C XP AND YP CONTAIN INFORMATION ABOUT THE CURVATURE OF THE                      
C CURVE AT THE GIVEN NODES,                                                     
C AND                                                                           
C S CONTAINS THE POLYGONAL ARCLENGTH OF THE CURVE.                              
      NM1 = N - 1                                                               
      NP1 = N + 1                                                               
C SET UP RIGHT HAND SIDES OF TRIDIAGONAL (WITH CORNER                           
C ELEMENTS) LINEAR SYSTEM FOR XP AND YP                                         
      DELX1 = X(2) - X(1)                                                       
      DELY1 = Y(2) - Y(1)                                                       
      DELS1 = SQRT(DELX1*DELX1+DELY1*DELY1)                                     
      DX1 = DELX1/DELS1                                                         
      DY1 = DELY1/DELS1                                                         
      XP(1) = DX1                                                               
      YP(1) = DY1                                                               
      TEMP(1) = DELS1                                                           
      S = DELS1                                                                 
      DO 10 I=2,N                                                               
        IP1 = I + 1                                                             
        IF (I.EQ.N) IP1 = 1                                                     
        DELX2 = X(IP1) - X(I)                                                   
        DELY2 = Y(IP1) - Y(I)                                                   
        DELS2 = SQRT(DELX2*DELX2+DELY2*DELY2)                                   
        DX2 = DELX2/DELS2                                                       
        DY2 = DELY2/DELS2                                                       
        XP(I) = DX2 - DX1                                                       
        YP(I) = DY2 - DY1                                                       
        TEMP(I) = DELS2                                                         
        DELX1 = DELX2                                                           
        DELY1 = DELY2                                                           
        DELS1 = DELS2                                                           
        DX1 = DX2                                                               
        DY1 = DY2                                                               
C ACCUMULATE POLYGONAL ARCLENGTH                                                
        S = S + DELS1                                                           
   10 CONTINUE                                                                  
      XP(1) = XP(1) - DX1                                                       
      YP(1) = YP(1) - DY1                                                       
C DENORMALIZE TENSION FACTOR                                                    
      SIGMAP = ABS(SIGMA)*FLOAT(N)/S                                            
C PERFORM FORWARD ELIMINATION ON TRIDIAGONAL SYSTEM                             
      DELS = SIGMAP*TEMP(N)                                                     
      EXPS = EXP(DELS)                                                          
      SINHS = .5*(EXPS-1./EXPS)                                                 
      SINHIN = 1./(TEMP(N)*SINHS)                                               
      DIAG1 = SINHIN*(DELS*.5*(EXPS+1./EXPS)-SINHS)                             
      DIAGIN = 1./DIAG1                                                         
      SPDIG1 = SINHIN*(SINHS-DELS)                                              
      SPDIAG = 0.                                                               
      DO 20 I=1,N                                                               
        DELS = SIGMAP*TEMP(I)                                                   
        EXPS = EXP(DELS)                                                        
        SINHS = .5*(EXPS-1./EXPS)                                               
        SINHIN = 1./(TEMP(I)*SINHS)                                             
        DIAG2 = SINHIN*(DELS*(.5*(EXPS+1./EXPS))-SINHS)                         
        IF (I.EQ.N) GO TO 30                                                    
        DIAGIN = 1./(DIAG1+DIAG2-SPDIAG*TEMP(I-1))                              
        XP(I) = DIAGIN*(XP(I)-SPDIAG*XP(I-1))                                   
        YP(I) = DIAGIN*(YP(I)-SPDIAG*YP(I-1))                                   
        TEMP(N+I) = -DIAGIN*TEMP(NM1+I)*SPDIAG                                  
        IF (I.EQ.1) TEMP(NP1) = -DIAGIN*SPDIG1                                  
        SPDIAG = SINHIN*(SINHS-DELS)                                            
        TEMP(I) = DIAGIN*SPDIAG                                                 
        DIAG1 = DIAG2                                                           
   20 CONTINUE                                                                  
   30 TEMP(NM1) = TEMP(N+NM1) - TEMP(NM1)                                       
      IF (N.EQ.2) GO TO 50                                                      
C PERFORM FIRST STEP OF BACK SUBSTITUTION                                       
      DO 40 I=3,N                                                               
        IBAK = NP1 - I                                                          
        XP(IBAK) = XP(IBAK) - TEMP(IBAK)*XP(IBAK+1)                             
        YP(IBAK) = YP(IBAK) - TEMP(IBAK)*YP(IBAK+1)                             
        TEMP(IBAK) = TEMP(N+IBAK) - TEMP(IBAK)*TEMP(IBAK+1)                     
   40 CONTINUE                                                                  
   50 XP(N) =                                                                   
     * (XP(N)-SPDIG1*XP(1)-SPDIAG*XP(NM1))/(DIAG1+DIAG2+SPDIG1*T                
     * EMP(1)+SPDIAG*TEMP(NM1))                                                 
      YP(N) =                                                                   
     * (YP(N)-SPDIG1*YP(1)-SPDIAG*YP(NM1))/(DIAG1+DIAG2+SPDIG1*T                
     * EMP(1)+SPDIAG*TEMP(NM1))                                                 
C PERFORM SECOND STEP OF BACK SUBSTITUTION                                      
      DO 60 I=1,NM1                                                             
        XP(I) = XP(I) + TEMP(I)*XP(N)                                           
        YP(I) = YP(I) + TEMP(I)*YP(N)                                           
   60 CONTINUE                                                                  
      RETURN                                                                    
      END                                                                       
      SUBROUTINE KURVP2(T, XS, YS, N, X, Y, XP, YP, S, SIGMA)                   
      INTEGER N                                                                 
      REAL T, XS, YS, X(N), Y(N), XP(N), YP(N), S, SIGMA                        
C THIS SUBROUTINE PERFORMS THE MAPPING OF POINTS IN THE                         
C INTERVAL (0.,1.) ONTO A CLOSED CURVE IN THE PLANE. THE                        
C SUBROUTINE KURVP1 SHOULD BE CALLED EARLIER TO DETERMINE                       
C CERTAIN NECESSARY PARAMETERS. THE RESULTING CURVE HAS A                       
C PARAMETRIC REPRESENTATIONBOTH OF WHOSE COMPONENTS ARE                         
C PERIODIC SPLINES UNDER TENSION AND FUNCTIONS OF THE POLY-                     
C GONAL ARCLENGTH PARAMETER.                                                    
C ON INPUT--                                                                    
C T CONTAINS A REAL VALUE OF ABSOLUTE VALUE LESS THAN OR                        
C EQUAL TO 1. TO BE MAPPED TO A POINT ON THE CURVE. THE                         
C SIGN OF T IS IGNORED AND THE INTERVAL (0.,1.) IS MAPPED                       
C ONTO THE ENTIRE CLOSED CURVE. IF T IS NEGATIVE THIS                           
C INDICATES THAT THE SUBROUTINE HAS BEEN CALLED PREVIOUSLY                      
C (WITH ALL OTHER INPUT VARIABLES UNALTERED) AND THAT                           
C THIS VALUE OF T EXCEEDS THE PREVIOUS VALUE IN ABSOLUTE                        
C VALUE. WITH SUCH INFORMATION THE SUBROUTINE IS ABLE TO                        
C MAP THE POINT MUCH MORE RAPIDLY. THUS IF THE USER SEEKS                       
C TO MAP A SEQUENCE OF POINTS ONTO THE SAME CURVE,                              
C EFFICIENCY IS GAINED BY ORDERING THE VALUES INCREASING                        
C IN MAGNITUDE AND SETTING THE SIGNS OF ALL BUT THE FIRST,                      
C NEGATIVE,                                                                     
C N CONTAINS THE NUMBER OF POINTS WHICH WERE INTERPOLATED                       
C TO DETERMINE THE CURVE,                                                       
C X AND Y ARE ARRAYS CONTAINING THE X- AND Y-COORDINATES                        
C OF THE INTERPOLATED POINTS,                                                   
C XP AND YP ARE THE ARRAYS OUTPUT FROM KURVP1 CONTAINING                        
C CURVATURE INFORMATION,                                                        
C S CONTAINS THE POLYGONAL ARCLENGTH OF THE CURVE,                              
C SIGMA CONTAINS THE TENSION FACTOR (ITS SIGN IS IGNORED).                      
C THE PARAMETERS N,X,Y,XP,YP,S AND SIGMA SHOULD BE INPUT                        
C UNALTERED FROM THE OUTPUT OF KURVP1.                                          
C ON OUTPUT--                                                                   
C XS AND YS CONTAIN THE X- AND Y-COORDINATES OF THE IMAGE                       
C POINT ON THE CURVE.                                                           
C T,N,X,Y,XP,YP,S AND SIGMA ARE UNALTERED.                                      
C DENORMALIZE SIGMA                                                             
      SIGMAP = ABS(SIGMA)*FLOAT(N)/S                                            
C STRETCH UNIT INTERVAL INTO ARCLENGTH DISTANCE                                 
      TN = ABS(T*S)                                                             
C FOR NEGATIVE T START SEARCH WHERE PREVIOUSLY TERMINATED,                      
C OTHERWISE START FROM BEGINNING                                                
      IF (T.LT.0.) GO TO 10                                                     
      I1 = 2                                                                    
      SUM = 0.                                                                  
   10 IF (I1.EQ.1) GO TO 50                                                     
C DETERMINE INTO WHICH SEGMENT TN IS MAPPED                                     
      DO 30 I=I1,N                                                              
        DELX = X(I) - X(I-1)                                                    
        DELY = Y(I) - Y(I-1)                                                    
        DELS = SQRT(DELX*DELX+DELY*DELY)                                        
        IF (SUM+DELS-TN) 20, 40, 40                                             
   20   SUM = SUM + DELS                                                        
   30 CONTINUE                                                                  
      I = 1                                                                     
      IM1 = N                                                                   
      DELS = S - SUM                                                            
      GO TO 50                                                                  
   40 IM1 = I - 1                                                               
C SET UP AND PERFORM INTERPOLATION                                              
   50 DEL1 = TN - SUM                                                           
      DEL2 = DELS - DEL1                                                        
      EXPS1 = EXP(SIGMAP*DEL1)                                                  
      SINHD1 = .5*(EXPS1-1./EXPS1)                                              
      EXPS = EXP(SIGMAP*DEL2)                                                   
      SINHD2 = .5*(EXPS-1./EXPS)                                                
      EXPS = EXPS1*EXPS                                                         
      SINHS = .5*(EXPS-1./EXPS)                                                 
      XS = (XP(I)*SINHD1+XP(IM1)*SINHD2)/SINHS +                                
     * ((X(I)-XP(I))*DEL1+(X(IM1)-XP(IM1))*DEL2)/DELS                           
      YS = (YP(I)*SINHD1+YP(IM1)*SINHD2)/SINHS +                                
     * ((Y(I)-YP(I))*DEL1+(Y(IM1)-YP(IM1))*DEL2)/DELS                           
      I1 = I                                                                    
      RETURN                                                                    
      END                                                                       
